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ABSTRACT 


A  software  package  has  been  assembled  (written  and  adapted  from  commercial  sources)  that  allows 
three  dimensional  atomic  level  modeling,  simulation  and  visualization  of  polycrystalline  samples 
subjected  to  various  stresses.  Sample  preparation,  is  accomplished  with  a  commercial  software 
package,  ADESH,  that  allows  initial  atomic  coordinates  to  be  prepared  conveniently  for  submission 
to  the  dynamics  portion  of  the  code.  The  code,  mdem.exe,  that  traces  the  atom  motions  was 
written  in  the  C  programming  language  for  convenient  portability  across  computing  platforms.  On 
an  Intel  486/66  class  machine,  it  can  provide  10  ps  of  history  for  a  system  containing  1896  atoms  in 
approximately  62  hours  of  computing  time.  Temperature  control  is  provided  by  a  thermostat  function 
known  to  preserve  the  statistical  integrity  of  the  dynamic  variables.  Several  techniques  have  been 
developed  for  visualization  of  the  output  data.  Preliminary  results  demonstrating  thermal  expansion, 
specific  heat,  self-diffusion  and  annealing  behavior  have  been  published. 


Task  Status 

As  the  attached  Man-loading  Summary  shows,  1544  of  the  1600  hours  contracted  for  Dr.  Helbig's 
time  under  Subcontract  CB0315  were  completed  by  the  end  of  the  contract  term.  Mr.  Bartelt 
completed  760  of  the  800  hours  that  had  been  proposed.  Thus,  of  the  anticipated  2400  man  hours, 
2304  or  96%  were  completed  and  billed.  The  difference  between  the  actual  and  anticipated  hours 
is  due  primarily  to  the  fact  that  Professor  Helbig  was  hired  as  a  full-time  employee  of  Rome  Lab 
approximately  90  calendar  days  before  the  scheduled  contract  completion  date.  With  the  submission 
of  this  report,  the  work  of  Subcontract  CB0315  is  concluded. 

Summary  of  Work  Completed 

Although  we  did  not  complete  all  the  work  we  had  hoped  to,  virtually  all  of  the  tasks  specified  in  the 
Statement  of  Work  were  completed.  The  chief  disappointment  resulted  from  an  overly  optimistic 
assessment  of  the  "electron  wind"  problem.  This  final  report  takes  as  its  six  Section  headings  the  six 
sentences  of  the  Statement  of  Work,  and  describes  how  each  of  the  specified  tasks  was  accomplished. 


1.  OVERVIEW 

"The  purpose  of  this  effort  is  to  develop  atomic  level  modeling  and  simulation 
software  that  characterizes  metallization  stripes  and  their  resistance  to  mass 
transport." 


From  the  outset,  it  was  clear  that  computing  power  will  not  be  available  in  the  near  term  to  simulate 
anything  approaching  a  complete  metallization  stripe.  The  best  current  fabrication  technology  cannot 


Page  2  of  19  -  Final  Report  for  CB03 1 5 


deliver  a  workable  stripe  with  a  volume  significantly  less  than  1  cubic  micron  containing  about  1012 
atoms.  The  best  current  computer  technology  can,  in  practical  times,  handle  computations  involving 
perhaps  106  atoms.  Wording  later  in  the  Statement  of  Work  (Section  2.)  makes  it  clear  that  subunits 
of  a  stripe,  consisting  of  a  few  crystalline  grains,  are  the  systems  of  interest  in  this  work.  In 
particular,  the  interfaces  or  grain  boundaries  between  two  such  crystallites  appear  to  be  the  most 
likely  routes  for  mass  transport  inside  a  metallic  system. 

At  the  conclusion  of  the  project,  a  software  package  has  been  assembled  (written  and  adapted  from 
commercial  sources)  that  allows  three  dimensional  atomic  level  modeling,  simulation  and  visualization 
of  polycrystalline  samples  subjected  to  various  stresses.  Sample  preparation  and  visualization  are 
handled  by  commercially  available  programs  as  described  in  Sections  2.  and  6.  below. 

The  program  mdem.exe,  compiled  from  its  C-language  source-code  modules,  was  written  to  follow 
efficiently  the  motions  of  assemblages  of  atoms  presented  to  it  as  a  file  listing  their  phase  space 
coordinates  (positions  and  velocities)  at  some  instant  and  species  identifiers.  The  program  reports 
its  results  in  the  form  of  files  that  record  a)  the  system  temperature  versus  time  and  b)  the  system 
phase  space  coordinates  versus  time.  The  output  files  containing  phase  space  coordinates  are  in  the 
same  format  as  the  input  file  making  it  convenient  to  chain  data  collection  runs  together. 

This  software  has  been  tested  and  debugged  on  a  several  model  systems.  Results  for  a  single  crystal 
sample  of  43  atoms  and  a  bi-crystal  sample  of  601  atoms  have  been  published1.  During  the  final 
weeks  of  the  contract,  a  bi-crystal  containing  1896  atoms  with  an  oblique  grain  boundary  was  studied. 
Both  the  601-atom  and  1896-atom  bi-crystal  samples  were  subjected  to  thermal  and  mechanical 
stresses  that  resulted  in  the  transport  of  atoms  to  locations  significantly  removed  from  their  original 
locations.  These  atomic  motions  were,  however,  minority  events;  the  vast  majority  of  atoms 
displayed  a  remarkable  tendency  to  retain  the  crystalline  structure  with  which  they  were  originally 
endowed. 

It  should  be  noted  that  for  both  of  these  bi-crystal  experiments,  no  external  constraints  were  imposed 
on  the  samples  that  would  have  artificially  maintained  the  integrity  of  their  crystal  structure. 
Nevertheless,  the  crystal  structure  was  retained  locally  (with  small  perturbations)  by  the  simple 
(Lennard-Jones)  force  model  used.  The  accumulation  of  small  perturbations  to  the  perfect  crystal 
structure  resulted  in  long  range  distortions  such  as  twisting  and  bending.  These  effects  are  evident 
from  a  comparison  of  Figures  1.  and  2.  which  show  orthographic  projections  of  the  initial 
configuration  of  the  601 -atom  sample  and  its  condition  after  it  had  stabilized  at  123  K.  The  atom 
layers  were  initially  arrayed  in  planes  parallel  with  the  x-y  plane.  The  grain  boundary  was  essentially 
coincident  with  the  plane  x  =  0.  The  evident  twisting  and  bending  are  due  to  a  large  initial  attractive 
force  between  the  two  crystalline  components,  to  the  misalignment  of  their  crystal  axes,  and  to  the 
asymmetry  of  the  two  crystallites;  that  on  the  left  contains  fewer  atoms  and  is  missing  most  of  its  z 
=  0  plane.  The  lone  atom  on  the  left  was  ejected  by  the  shock  wave  produced  by  the  collision  of  the 
two  crystallites  and  then  captured  as  an  adsorbed  particle. 
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2.  SAMPLE  PREPARATION 


", Acceptable  software  shall  include  tools  for  constructing  thin  film  metallization 
consisting  of  a  minimum  of  three  grains  that  accurately  simulates  microelectronic 
metallization  microstructure  in  three  dimensions.  ” 


As  indicated  in  Section  1.,  the  strategy  that  evolved  to  create  the  simulation  package  called  for  in  the 
Statement  of  Work  was  to  separate  the  three  constituent  tasks;  1)  constructing  suitable  arrays  of 
initial  atomic  coordinates,  2)  following  the  motions  of  the  atoms  as  they  moved  from  their  initial 
positions,  and  3)  analyzing  the  results. 

The  first  task  amounts  to  sample  preparation,  and  is  chiefly  a  geometric  one;  modeling  in  the  sense 
of  molding  or  sculpting  the  atoms  into  sensible  arrangements  to  simulate  portions  of  a  polycrystalline 
solid.  Two  complete  solutions  to  the  sample  preparation  problem  were  found  in  the  course  of  the 
work.  A  working  sketch  called  MAKECRYS  .  BAS  was  constructed  in  QuickBasic,  and  with  it  the  601 
atom  bicrystal  was  prepared  for  dynamic  study.  Much  of  this  sketch  was  reduced  to  C-code,  but 
development  was  deferred  while  the  second  task,  that  of  the  particle  dynamics,  was  solved.  Our 
expectations  that  this  would  provide  us  with  a  better  understanding  of  the  overall  problem  were 
fulfilled.  Before  we  returned  to  the  sample  preparation  problem,  a  commercial  software  package, 
(ADESH©  from  CASA)  was  identified  that  allows  initial  atomic  coordinates  to  be  prepared 
conveniently  for  submission  to  the  dynamics  portion  of  the  code. 

The  1896-atom  bi-crystal  mentioned  in  Section  1.  was  created  by  ADESH.  It  is  shown  in  near- 
orthographic  projection  in  Figure  3.  Figure  4.  shows  the  same  system  after  annealing  at  10  K. 
Further  discussion  of  this  system  appears  in  Section  4. 

It  is  not  difficult  to  generate  the  coordinates  of  atoms  that  represent  a  crystallite  of  desired  lattice 
constant  and  orientation.  One  merely  defines  a  suitably  positioned  and  oriented  set  of  basis  vectors, 
and  then  systematically  multiplies  them  by  integers.  If  the  crystal  requires  more  than  one  atom  per 
unit  cell,  that  presents  only  small  complication;  a  set  of  basis  vectors  internal  to  the  unit  cell  must  be 
defined,  and  used  to  add  the  appropriate  additional  atoms  to  each  unit  cell. 

The  chief  difficulty  lies  in  terminating  the  list  of  atomic  coordinates  so  that  they  lie  within  some 
desired  geometric  bounds.  In  the  QuickBasic  sketch  mentioned  above,  we  provided  for  crystallites 
bounded  by  convex,  but  otherwise  arbitrary  polyhedrons  with  up  to  twenty  faces,  each  face 
containing  up  to  ten  edges.  Within  each  polyhedron  the  user  could  specify  the  coordinates  of  one 
atom,  and  a  set  of  basis  vectors  to  fix  the  crystal  orientation  and  lattice  constant.  The  program  then 
filled  the  polyhedron  by  systematically  adding  atoms  coordinates  and  checking  for  each  new  trial 
coordinate  set  that  it  was  inside  the  bounding  polyhedron.  The  search  for  additional  coordinates  was 
continued  out  to  a  rectangular  parallelepiped,  defined  by  the  basis  vectors,  that  circumscribed  the 
desired  polyhedron. 
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The  ADESH  code  presumably  goes  through  a  similar  process  to  generate  sets  of  atomic  coordinates, 
but  the  authors  simplified  their  task  by  restricting  the  bounding  solids  to  rectangular  parallelepipeds, 
or  segments  (including  the  full  figure)  of  a  cylinder  or  sphere.  Although  this  seemed  at  first  to  be 
restrictive,  we  were  able  to  construct  a  bicrystal  sample  that  included  a  grain  boundary  oblique  to  the 
bounding  surfaces  as  Figure  3.  illustrates.  We  have  not,  in  fact,  found  ways  to  construct  crystallites 
using  ADESH  that  would  combine  to  form  a  truly  three  dimensional  triple  point  joining  three 
crystallites.  If  CASA,  Inc.  or  another  commercial  vendor  does  not  provide  a  solution  to  this  problem, 
it  will  be  necessary  to  complete  the  development  of  MAKECRYS .  BAS  into  C-code. 


3.  PHYSICAL  CONSTANTS 

"The  software  shall  be  capable  of  determining  accurate  physical  constants  for  the 
material  being  modeled,  such  as  Young's  modulus,  thermal  coefficient  of 
expansion,  thermal  conductivity,  diffusion  coefficients,  etc  " 


Our  treatment  of  physical  constants  proceeded  in  this  way.  The  two  free  parameters  in  the  Lennard- 
Jones  potential  function,  which  was  used  to  model  the  atomic  interactions,  were  chosen  to  reproduce 
the  known  lattice  parameter  and  Young's  modulus  for  aluminum. 

The  resulting  thermal  coefficient  of  expansion,  a,  was  then  determined  for  a  43-atom  spherical  cluster 
from  the  slope  of  the  fractional  change  of  its  average  maximum  dimension  versus  its  temperature. 
The  value  of  a  obtained  in  this  way  was  larger  than  measured  values  for  aluminum  by  a  factor  of 
about  20.  This  was  attributed  to  the  fact  that  the  Lennard-Jones  binding  energy  was  modeled  on  the 
engineering  value  of  Young's  modulus  for  Al,  not  on  the  larger  value  that  would  pertain  to  an  ideal 
crystal.  Had  the  larger  value  been  used,  a  would  have  been  in  closer  accord  with  measurement.  It 
is  pertinent  that  the  ADESH  program  mentioned  above  uses  a  Lennard-Jones  binding  energy  for  Al 
that  is  eight  times  larger  than  that  used  in  our  work.  Their  value  is  in  accord  with  the  heat  of 
formation  of  aluminum. 

At  best,  however,  the  Lennard-Jones  potential  has  not  proven  very  effective  as  a  model  of  metallic 
bonding,  and  no  significant  time  was  spent  exploring  the  details  of  this  potential  model.  The 
simulation  code  will  easily  accept  alternative  potential  functions;  we  are  looking  forward  to 
implementing  more  sophisticated  interaction  models. 

The  self-diffusion  coefficient  for  Al  was  also  determined.  This  was  done  for  the  601-atom  bi-crystal 
by  recording  the  mean-square  displacement  of  its  atoms  after  it  had  reached  a  stable  mean  kinetic 
energy  corresponding  to  a  temperature  of  123  K.  Not  surprisingly  for  this  temperature,  the  self- 
diffusion  coefficient  was  found,  within  the  error  imposed  by  thermal  fluctuations  for  such  a  small 
system,  to  be  zero.  This  result  is  consistent  with  the  stability  of  the  system  that  was  observed  using 
a  variety  of  visualization  techniques.  Experiments  at  higher  temperatures  and  with  the  presence  of 
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vacancies  have  not  yet  been  made,  but  the  capability  of  doing  so  has  been  demonstrated. 

A  third  physical  constant  was  determined,  namely  the  specific  heat  of  the  43 -atom  cluster  mentioned 
above.  The  expected  value  for  a  system  constrained  by  the  laws  of  classical  mechanics  is  that  of 
Dulong  and  Petit,  and,  within  the  uncertainty  of  thermal  fluctuations,  that  was  the  value  found. 

No  measurements  have  yet  been  made  of  thermal  conductivity,  but  the  thermostat  function  contained 
in  the  simulation  could  easily  be  set  to  different  values  at  opposite  ends  of  a  sample  and  monitored 
for  the  time  averaged  injection  (at  the  hot  end)  and  extraction  (at  the  cold  end)  of  energy.  These, 
together  with  the  sample  dimensions  provide  the  necessary  data  to  compute  the  thermal  conductivity. 


4.  MECHANICAL,  THERMAL,  ELECTRICAL  STRESSES 

"  The  software  shall  provide  introducing  stress,  temperature,  and  thermal  and 
electrical  gradients  into  the  simulated  metallization  stripe  " 


Mechanical  stress  can  be  introduced  into  the  simulated  atomic  system  either  statically  or  dynamically. 
Static  stress  may  be  introduced  by  fixing  the  positions  of  atoms  at  opposite  ends  of  the  sample  in  such 
a  way  that  they  cannot  all  be  part  of  a  perfect  crystal.  If,  for  example,  the  entire  collection  of  surface 
atoms  is  fixed,  the  interior  regions  will  be  subjected  to  either  positive  or  negative  pressures  depending 
on  the  temperature  of  the  system.  Similarly,  atoms  on  opposite  faces  can  be  fixed  in  such  a  way  that 
sheer,  torsional,  compressive  or  tensile  stresses  (or  combinations  of  them)  are  imposed.  Although 
this  facility  is  present,  no  experiments  have  yet  been  conducted  employing  static  mechanical  stresses. 

Dynamic  mechanical  stress  is  necessarily  imposed  when  a  system  is  prepared  in  a  non-equilibrium 
configuration,  and  then  allowed  to  relax  via  molecular  dynamics  simulation.  Since  it  is  virtually 
impossible  to  select  the  initial  phase  space  coordinates  for  a  system  of  any  complexity  so  that  it  will 
be  in  an  equilibrium  configuration,  such  dynamic  stresses  are  unavoidable.  They  can  be  minimized, 
for  example,  by  utilizing  a  non-dynamic  technique  to  relax  the  system  before  beginning  molecular 
dynamics  simulation.  The  ADESH  program  uses  Monte  Carlo  methods  to  do  just  this. 

Alternatively,  one  can  rely  on  a  thermostat  function  to  cushion  the  effects  of  initial,  non-equilibrium 
conditions.  This  is  successful  provided  the  initial  conditions  are  reasonably  close  to  equilibrium.  The 
warped  condition  of  the  601 -atom  sample  pictured  in  Figure  2  and  its  having  annealed  from  a 
bicrystal  with  a  grain  boundary  to  an  imperfect  single  crystal  is  due  in  large  part  to  dynamical  stresses 
that  developed  as  the  two  crystallites  collided.  The  temperature  record  of  this  encounter  is  shown 
in  Figure  5.  When  the  same  system  was  allowed  to  relax  with  a  thermostat  function  active  and  set 
at  1  K,  the  peak  temperature  reached  was  3.5  K,  and  thereafter  the  system  temperature  settled  down 
within  10  ps  to  acceptably  small  fluctuations  about  1  K  as  shown  in  Figure  6.  Importantly,  the  grain 
boundary  was  preserved,  so  that  this  system  can  now  serve  as  the  starting  point  for  studies  of  mass 
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transport  along  (and  perpendicular  to)  a  grain  boundary  with  at  least  some  thermodynamic  stability. 
Whether  that  stability  will  survive  heating  to  interestingly  high  temperatures,  let  alone 
electromigration  stresses,  remains  to  be  seen. 

Several  experiments  of  this  nature  were  conducted,  and  it  was  found  that  the  thermostat  function  can 
become  unstable  if  one  attempts  to  change  the  system  temperature  too  rapidly.  With  2  fs  integration 
time  steps,  temperature  changes  as  large  as  10  K  are  successful.  Figure  7.  shows  the  well-behaved 
temperature  evolution  of  the  1896-atom  bicrystal  as  it  evolved  from  rest  (0  K)  as  shown  in  Figure  3. 
for  a  time  interval  of  10  ps.  Figure  4.  shows  its  final  configuration  at  10  K.  For  larger  temperature 
changes,  the  system  mean  kinetic  energy  will  fluctuate  with  large  and  often  catastrophic  amplitude. 
Figure  8.  shows  the  result  of  attempting  to  raise  the  temperature  of  the  1896-atom  sample  from  10 
K  to  100  K.  Clearly,  the  thermostat  was  unable  to  cope  with  this  large  temperature  change  using  the 
same  integration  time-step  and  time  constant  as  were  used  to  obtain  the  result  shown  in  Figure  7. 

From  the  foregoing  discussion,  it  is  apparent  that  temperature  control  is  provided  for  in  the  simulation 
code.  We  have  implemented  the  Nose-Hoover  thermostat  function  as  described  by  Valkealahti  and 
Manninen3  to  accomplish  this.  Although  no  experiments  have  yet  been  done  with  thermal  gradients, 
by  setting  the  temperatures  of  different  groups  of  atoms  to  different  temperatures,  such  gradients  can 
easily  be  created  and  studied  as  mentioned  in  Section  3  in  the  discussion  of  thermal  conductivity. 

Electric  fields  (i.e.,  voltage  gradients)  can  easily  be  introduced  as  additions  to  the  force  function  that 
operates  on  each  atom  in  the  function  move_atoms  in  the  C-code  module  mdemutil .  c.  A 
uniform  electric  field  acting  on  an  atom  (considered  as  a  positive  ion  due  to  its  sharing  of  non- 
localized  electrons  with  the  crystal)  would  be  represented  as  a  constant  additional  force  parallel  with 
the  direction  of  the  electric  field.  A  variation  of  this  notion  is,  in  fact,  a  standard  procedure  in 
continuum  models  of  electromigration  where  the  lattice  ions  are  endowed  with  an  effective  charge, 
Z*,  taken  to  negative  in  order  that  the  atoms  should  move  (in  accordance  with  observation)  anti¬ 
parallel  with  the  electric  field  (in  the  direction  of  the  electron  motion).  In  our  view,  this  ploy  to 
attribute  the  ion  motion  directly  to  the  electric  field,  rather  than  to  the  action  of  electrons  as  the 
agents  of  the  electric  field,  is  unwarranted  in  a  model  intended  to  reflect  as  accurately  as  possible  all 
of  the  influences  at  work  at  the  atomic  scale. 

Unfortunately,  we  have  not  yet  been  able  to  introduce  a  satisfactory  model  of  electron-ion  interaction 
(energy  and  momentum  transfer)  into  the  simulation.  Indeed,  to  our  knowledge,  no  one  has  achieved 
this.  As  mentioned  in  our  4th  quarterly  report,  a  standard  reference  in  the  field,  Ziman's  "Electrom 
and  Phonons",  declares 

"Let  it  be  said,  at  once,  that  there  is  no  rigorous  theory  of  the  transport  coefficients. 

The  relationship  between  the  electric  current  and  the  electric  field  in  a  metal  at  a  fixed 
temperature  has  never  been  formally  deduced  from  the  equations  of  motion  of  the 
electrons,  classical  or  quantal." 

We  are,  nonetheless,  trying  to  construct  an  electron-ion  interaction  model  that  will  provide  a 
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consistent  picture  of  momentum  transfer  from  the  electrons  to  the  ions  (the  "electron  wind")  and  of 
excess  energy  transfer  from  the  electrons  to  the  ions  (I2R,  or  Joule  heating).  Our  current  thinking  on 
this  subject  centers  on  the  following  notions.  Electrons  and  ions  in  a  conductor  are  constantly 
exchanging  energy  as  a  result  of  collision  processes.  At  thermal  equilibrium,  and  in  the  absence  of 
an  electric  current,  the  transfer  of  energy  from  the  electrons  to  the  lattice  and  to  the  electrons  from 
the  lattice  must  proceed  at  nearly  identical  rates.  The  difference  appears  to  be  I2R  heating  of  the 
lattice  by  thermally  driven  current  fluctuations  (Johnson  noise).  One  might  suppose  that  this  excess 
energy  delivered  to  the  lattice  escapes  from  the  conductor  radiatively.  For  this  picture  to  be 
consistent,  the  cooling  of  the  electrons  due  to  Johnson  noise  must  be  offset  by  the  energy  absorption 
processes  that  maintain  the  conductor  in  thermal  equilibrium  with  its  environment,  e.g!,  the  absorption 
of  radiation  from  that  environment.  These  are  the  background  activities  upon  which  is  superposed 
the  tiny  electron  drift  that  constitutes  even  a  large  DC  current  (i.e.,  the  electron  wind).  A  model  of 
electron-ion  interactions  that  accounts  for  these  processes  is  what  we  are  presently  seeking. 


5.  PROGRAMMING  LANGUAGE 

"  The  source  code  shall  be  delivered  to  Rome  Laboratory  for  future  development, 
shall  be  written  in  the  AN SII  Standard  Version  of  the  C programming  language, 
and  shall  build  upon  software  previously  developed  at  Rome  Laboratory  " 


The  simulation  software  generated  under  this  contract  is  an  outgrowth  of  software  developed  at 
Rome  Lab  during  the  Summer  of  1993.  It  has  been  delivered  to  Rome  Laboratory  and  is  under 
continuing  development  in  the  very  real  sense  that  the  principal  investigator  is  now  a  full-time 
employee  of  the  Lab  and  is  devoting  most  of  his  effort  to  improving  and  extending  the  software. 

The  executable  program,  mdem.exe,  that  follows  the  motions  of  the  atoms  comprising  the  sample 
system  was  compiled  from  four  modules  written  in  the  ANSI!  Standard  C  programming  language. 
The  latest  version  of  the  code  was  compiled  as  32-bit  DOS  code  to  run  on  a  486-class  personal 
computer.  The  compiler  chosen  for  this  work  was  the  WATCOM  C/C'"'  compiler  (version  10).  It 
has  the  capability  to  compile  and  optimize  the  same  source  code  for  other  operating  systems  including 
OS-2  and  Windows  NT.  Compilation  for  UNIX  based  machines  should  pose  no  problem  because 
the  WATCOM  integrated  development  environment  will  screen  the  source  code  for  ANSI! 
standardization. 

The  four  modules  are  mdem.c,  mdemio.c,  mdemutil.c,  and  mdemmem.c.  They  are 
accompanied  by  so-called  "header  files"  of  the  same  name  but  with  the  extension  .  h  that  contain 
declarations  of  external  variables  and  functions  from  the  standard  library. 

The  module  mdem.c  contains  the  function  named  "main"  that  the  C-compiler  treats  as  the  control 
module.  It  accepts  command-line  arguments  that  determine  the  data  file  in  which  the  phase  space 
coordinates  are  to  be  found,  the  number  of  integration  steps  to  be  made  (nsteps),  the  integration 
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time-step,  the  thermostat  function  initial  value,  the  desired  system  temperature  and  the  number  of 
integration  steps  between  the  output  of  phase  space  files.  Its  chief  function  is  to  execute  a  loop  in 
which  it  calls  the  function  move_atoms  nsteps  times.  It  also  keeps  track  of  elapsed  time,  moves 
array  pointers,  integrates  the  thermostat  function,  check  for  user  interrupts,  writes  useful  status 
information  to  the  screen  and  writes  data  files  to  disk  within  this  loop.  Prior  to  entering  the  loop,  it 
calls  functions  from  the  other  three  modules  to  perform  various  initializations,  and  after  the  loop 
terminates,  it  closes  the  program  in  an  orderly  fashion. 

The  module  mdemio .  c  contains  input/output  functions  used  by  mdem.  c  to  read  various  sorts  of 
data  from  the  initial  data  file,  to  acquire  date/time  information  from  the  operating  system  (for  use  in 
naming  output  files)  and  to  write  the  phase  space  coordinates  to  disk. 

The  important  content  of  the  module  mdemutil .  c  is  the  function  mo  ve_a  toms.  It  performs  the 
numerical  integrations  whereby  each  atom  is  moved  in  accord  with  its  current  phase  space 
coordinates  (position  and  momentum),  the  forces  exerted  on  it  by  all  the  other  atoms  within  its  sphere 
of  influence  and  the  force  that  the  thermostat  function  imposes  to  control  the  temperature.  This 
module  also  contains  functions  to  compute  the  system  kinetic  and  potential  energies,  and  some 
incidental  functions  dealing  with  the  system  linear  and  angular  momentum  (that  are  presently  unused). 

The  module  mdemmem.  c  contains  functions  that  allocate  and  de-allocate  memory  for  various  scalar, 
vector  and  tensor  quantities.  Some  of  these  are  currently  unused.  The  data  structure  of  primary 
importance  is  the  collection  of  phase  space  coordinates  and  a  parallel  array  containing  the  tags  that 
identify  the  species  of  each  atom.  Currently,  only  collections  of  identical  atoms  have  been  studied, 
but  the  tag  array  provides  for  differing  interactions  among  several  species. 


6.  VISUALIZATION 

"A  graphical  interface  for  constructing  the  model  and  interpreting  the  output  shall 
also  be  considered.” 


Several  means  were  explored  to  provide  visualization  of  the  data  that  were  collected  from  the 
simulated  experiments  described  in  the  previous  sections.  Figures  1.  through  4.  exemplify  the  most 
satisfactory  two-dimensional  graphical  presentations  of  the  atomic  configurations  studied.  They  were 
prepared  using  the  commercial  data  analysis  program  Origin  from  Microcal,  Inc.  Origin  is  also 
capable  of  displaying  atomic  configurations  as  two-dimensional  projections  (true  perspective  or 
orthographic)  of  their  three-dimensional  positions.  As  a  pertinent  example.  Figure  9.  is  a  perspective 
projection  of  the  1896-atom  sample  in  its  final,  10  K  condition. 

In  addition  to  Origin,  code  was  written  in  QuickBasic  to  produce  three-dimensional  perspective  views 
of  the  samples,  and  the  commercial  package  SigmaPlot  from  Jandel,  Inc.  was  evaluated.  The 
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drawback  with  QuickBasic  is  that  it  is  restricted  to  VGA  (640  X  480  pixels)  graphics.  This  is 
unacceptable  compared  with  the  1024  X  768  graphics  now  generally  available  even  on  PC-class 
computers.  SigmaPlot  was  rejected  because  it  cannot  render  a  three-dimensional  view  of  two  sets 
of  differently  colored  (or  shaded)  objects  that  properly  preserves  the  "back-to-ffont"  order  of  all 
objects.  The  result  is  that  objects  which  should  be  hidden  by  others  closer  to  the  viewpoint,  are  not. 


Finally,  a  cinematic  presentation  was  prepared  of  the  evolution  of  the  601 -atom  system  from  its  initial 
condition  (Figure  1.)  to  its  final  condition  (Figure  2.)  by  placing  95  intermediate  position  plots  (like 
those  of  Figures  1.  and  2.)  successively  on  the  video  screen.  In  the  interests  of  a  smooth  visual  effect, 
the  previous  plots  were  not  erased,  so  the  effect  was  to  see  the  paths  of  the  601  atoms  develop 
simultaneously.  By  watching  this  video  presentation  closely  and  repetitively,  it  was  possible  to  detect 
the  three  surges  of  kinetic  energy  that  are  visible  in  the  temperature  record  of  Figure  5;  the  initial 
impact  of  the  two  crystallites  is  followed  by  three  cycles  of  a  bending-mode  oscillation  before 
thermalization  sets  in. 

This  code  was  written  in  QuickBasic.  The  VGA  graphics  limitation  was  tolerable  for  this  application 
which  represented  each  atom  position  as  a  single  pixel.  The  projected  atom  positions  were  reduced 
to  integer  representation  to  allow  maximum  speed  of  data  handling  and  to  maximize  the  number  of 
points  that  could  be  held  in  RAM. 
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Figure  1  Simulated  Al  bi-crystal  at  rest  (0  K). 
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Figure  2.  Simulated  Al  crystallite,  after  12  ps  (123  K). 
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X  Axis  (Angstroms) 

Figure  3.  1896-atom  bi-crystal  from  ADESH  (Initial  Condition) 
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X  Axis 

Figure  4.  1896-atom  bi-crystal  (10  K) 
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Figure  5  Temperature  vs.  Time  for  601  atom 
bi-crystal.  Auxiliary  lines  show  mean  temperature, 
rms  fluctuation  and  stat.  mech.  expected  fluctuation. 
Total  time  ~  13  ps. 
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Figure  6  Temperature  versus  time 

for  the  601 -atom  system  with  the  thermostat 
function  active  and  set  to  1  K. 
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Figure  7  Temperature  versus  time  for  1896-atom 
bicrystal  starting  from  0  K  and  with 
thermostat  set  at  10  K. 
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Figure  9.  1896-atom  bi-crystal  (10  K) 


MISSION 

OF 

ROME  LABORATORY 

Mission.  The  mission  of  Rome  Laboratory  is  to  advance  the  science  and 
technologies  of  command,  control,  communications  and  intelligence  and  to 
transition  them  into  systems  to  meet  customer  needs.  To  achieve  this, 
Rome  Lab: 

a.  Conducts  vigorous  research,  development  and  test  programs  in  all 
applicable  technologies; 

b.  Transitions  technology  to  current  and  future  systems  to  improve 
operational  capability,  readiness,  and  supportability; 

c.  Provides  a  full  range  of  technical  support  to  Air  Force  Materiel 
Command  product  centers  and  other  Air  Force  organizations; 

d.  Promotes  transfer  of  technology  to  the  private  sector; 

e.  Maintains  leading  edge  technological  expertise  in  the  areas  of 
surveillance,  communications,  command  and  control,  intelligence,  reliability 
science,  electro-magnetic  technology,  photonics,  signal  processing,  and 
computational  science. 

The  thrust  areas  of  technical  competence  include:  Surveillance, 
Communications,  Command  and  Control,  Intelligence,  Signal  Processing, 
Computer  Science  and  Technology,  Electromagnetic  Technology, 
Photonics  and  Reliability  Sciences. 


